(12) INTERNATIONAL APPLICATION PUBLISHED UNDER THE PATENT COOPERATION TREATY (PCT) 



(19) World Intellectual Property Organization 
International Bureau 



lillliiillllllillllllllilllillllllliillira 



(43) International Publication Date 
22 May 2003 (22.05.2003) 



PCT 



(10) International Publication Number 

WO 03/041584 A2 



(51) International Patent Classification^: A61B 6/00 
(21) International Application Number: PCT/US02/3604S 



(22) International Filing Date: 

12 November 2002 (12.11.2002) 



(25) Filing Language: 

(26) Publication Language: 



English 
Englisli 



(30) Priority Data; 

10A)10,773 13 Navember200l (13.11.2001) US 

(71) Applicants: KONINKLIJKE PHILIPS ELECTRON- 
ICS NV [NL/NL]; Gracnewoudseweg 1, NL-5621 BA 
J Eindhoven (NL). PHILIPS MEDICAL SYSTEMS 
: (CLEVELAND), INC. [US/US]; 595 Miner road, Cleve- 
; land, OH 44143 (US). 



(74) Agents: LUNDIN, Thomas, M,; Philips Medical Systems 
(Cleveland), Inc., 595 Miner Road, Cleveland, OH 44143 
etal. (US). 

(81) Designated State (national): JP. 

(84) Designated States (regional): European patent (AT, BE, 
BG, CH, CY, CZ, DE, DK, EE, ES, FI, FR, GB, GR, IE. IT, 
LU, MC, NI-. PT, SE, SK, TR). 

Published: 

— : without international search report and to be republished 
upon receipt of thai report 

For two-letter codes and other abbreviations, refer to the "Guid- 
ance Notes on Codes and Abbreviations " appearing at the begin- 
ning of each regular issue of the PCT Gazette. 



■ (72) Inventors: SURl, Jasjit, S.; 655 Dewitt Drive, Highland 

s Heights, OH 44143 (US). LIU, Kecheng; 32701 S. Round- 

: head Drive, Solon, OH 44139 (US). WU, Dee, H.; 2931 

; Huntington Road, Shaker Heights, OH 44120 (US). 



= (54) Title: ANGIOGRAPHY MmnOD AND APPARATUS 




^ (57) Abstract: A two-dimensional slice formed of pixels (37(5) is extracted fiiom the angiographic image (76) after enhancing the 
^ vessel edges by second order spatial differentialion (78). Imaged vascular slructures in the slice are Oood-lillod (384). The edges of 
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ANGIOGRAPIT? METHOD AND APPARATUS 

Background o£ the Invention 

The present invention relates to the medical imaging 
5 arts. It particularly relates to angiography using the 
magnetic resonance imaging (MRI) and computed tomography (CT) 
medical imaging techniques, and will be described with 
particular reference thereto. However, the invention will also 
find application in conjunction with other three-dimensional 
10 imaging modalities as well as in other imaging arts in which 
thin structures or networks with overlapping or furcated 
portions are advantageously differentiated from extraneous 
imaged structures and background noise or tracked in three 
dimensions . 

15 Catastrophic medical events such as heart attacks and 

strokes that result from underlying vascular problems are a 
leading cause of death in the United States. Plaque buildup on 
the inside of the vascular walls can lead to strokes, coronary 
heart disease, and other medical conditions such as vessel 

20 rupture. Many Americans suffer from chronic vascular diseases 
which degrade quality of life. 

Angiography relates to the imaging of blood vessels 
and blood vessel systems, and as such is a powerful medical 
diagnostic for identifying and tracking vascular diseases. 

25 Angiography enables improved surgical planning and treatment, 
improved diagnosis and convenient non- invasive monitoring of 
chronic vascular diseases, and can provide an early indication 
of potentially fatal conditions such as aneurysms and blood 
clots. The ability of certain types of angiography to 

3.0 accurately characterize the vessel lumen is particularly 
valuable for diagnosing plaque buildup on the vessel walls. 

Angiography is performed using a number of different 
medical imaging modalities, including biplane X-ray/DSA, 
magnetic resonance (MR) , computed tomography (CT) , ultrasound, 

35 and various combinations of these techniques. Two-dimensional 
or three-dimensional angiographic data can be acquired 
depending upon the medical imaging modality and the selected 
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densities such as in the brain. Bifurcation points, tortuous 
or occluded vessels, vessel overlaps, intertwined vessels, 
partial volume averaging and other imaging artifacts, and 
vessel gaps can act alone or in various combinations to produce 
5 localized disjoints of the vascular path (or localized 
disjoints of the angiographic image of the vascular path) which 
prevent successful vessel tracking. Furthermore, at vessel 
overlaps the wrong vascular system may be tracked. For 
example, a tracking system following the arterial branch A of 
10 FIGURE 1 could fail and begin tracking the venous branch V at 
any of the crossing points AV, VA. 

The present invention contemplates an improved 
angiographic method and apparatus which overcomes the 
aforementioned limitations and others. 

15 

Summary of the Invention 
According to one aspect of the invention, a method is 
provided for characterizing a vascular system in a three- 
dimensional angiographic image comprised of voxels. A two- 

20 dimensional slice formed of pixels is extracted from the 
angiographic image. Vessel center points are identified in the 
two-dimensional slice. The extracting and identifying are 
repeated for a plurality of slices to generate a plurality of 
vessel center points that are representative of the vascular 

25 system. Segmenting, tracking, extracting, enhancing, or 
identifying of vascular information contained in the 
angiographic image is performed using the vessel center points 
as operative inputs. 

According to another aspect of the invention, an 

30 apparatus is disclosed for characterizing a vascular system. An 
imaging means is provided for acquiring imaging data from at 
least a portion of a patient. A reconstruction means is 
provided for reconstructing the acquired imaging data to 
generate a three-dimensional image representation comprised of 

35 voxels and including enhanced vascular contrast. A means is 
provided for extracting from the image representation a 
two-dimensional slice formed of pixels. A means is provided for 

- 4 - 
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identifying vessel center points in the two-dimensional slice, 
A means is provided for iteratively invoking the extracting and 
identifying means for a plurality of slices to generate a 
plurality of vessel center points that are representative of 
5 the vascular system. A means is provided for segmenting, 
tracking, extracting, enhancing, or identifying vascular 
information contained in the angiographic image using the 
vessel center points as operative inputs. 

One advantage of the present invention is that it is 
10 a global technique V7hich overcomes segmentation failures often 
encountered at vessel disjoints and overlaps by localized 
techniques . 

Another advantage of the present invention is that it 
is compatible with both white blood angiography and black blood 
15 angiography. 

Another advantage of the present invention is that it 
provides rapid and accurate vessel boundary estimation using 
propagation of geometric contours. The propagating can 
advantageously incorporate the gray scale image information 
20 through fuzzy membership classification of image elements in 
the neighborhood of the contour. 

Another advantage of the present invention is that it 
provides locations of furcation and vessel overlap points 
globally throughout the angiographic volume. This information 
25 can be used by vessel trackers or other vessel segmentation 
systems to improve accuracy and speed. 

Yet another advantage of the present invention is 
that it improves tracking speed and accuracy by providing a 
plurality of vessel center points or tags that are 
30 representative of the global vascular system. 

Still yet another advantage of the present invention 
is that it advantageously retains the sharp vascular edges and 
accurate vessel lumen information typically achieved by black 
blood data acquisition during the vessel segmentation 
35 processing. 

Numerous additional advantages and benefits of the 
present invention will become apparent to those of ordinary 
skill in the art upon reading the following detailed 

- 5 - 
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description of the preferred embodiment. 

Brief Description of the Drawinga 

The invention may take form in various components and 
5 arrangements of components, and in various steps and 
arrangements of steps. The drawings are only for the purpose 
of illustrating preferred embodiments and are not to be 
construed as limiting the invention. 

FIGURE 1 shows a schematic drawing of a pair of 
10 intertwined vascular structures that have several bifurcation 
points; 

FIGUKE 2 shows a schematic drawing of a vessel 
crossing with typical angiographic slice images superimposed, - 

FIGURE 3 A shows a typical angiographic slice in the 
15 overlap region of FIGURE 2 with the vessels artificially 
differentiated; 

FIGURE 3B shows a typical angiographic slice in the 
overlap region of FIGURE 2 without the artificial 
dif f erentiation; 
20 FIGURE 4 shows an exemplary magnetic resonance 
angiography (MRA) system formed in accordance with an 
embodiment of the invention; 

FIGURE 5 shows an overview of an exemplary 
angiographic segmentation method fqrmed in accordance with an 
25 embodiment of the invention; 

FIGURE 6A schematically shows a plurality of vessel 
centers that are representative of a bifurcated vascular 
branch; 

FIGURE 6B schematically shows a plurality of 
30 estimated vessel boundaries corresponding to the vessel centers 
of FIGURE 6A; 

FIGURE 6C schematically shows a representation of the 
bifurcated vascular branch reconstructed using the vessel 
centers and boundaries of FIGURES 6A and 6B; 
35 FIGURE 7 shows an exemplary embodiment of the pre- 

processor of the segmenter of FIGURE 5; 

FIGURE 8 shows an exemplary embodiment of the pseudo- 
white blood angiographic (pseudo-WBA) processor of FIGURE 7; 

- 6 - 
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FIGURE 9 shows an exemplary embodiment of the bone- 
air-tnuscle mask processor of FIGURE 8; 

FIGURE 10 shows an exemplary embodiment of the class 
assignment processor of FIGURE 9; 
5 FIGURE 11 shows an exemplary embodiment of the three- 

dimensional high curvature feature processor of FIGURE 5; 

FIGURE 12 shows an exemplary embodiment of the 
decomposition of the Gaussian derivative kernel into 
cosinusoidal and sinusoidal components; 
10 FIGURE 13 shows an exemplary time-of -flight white 

blood angiographic image of the carotid area that includes five 
isolated vessels as well as one pair of vessels emerging from 
a bifurcation; 

FIGURE 14 shows the edge volume corresponding to the 
15 image of FIGURE 13; 

FIGURES 15 and 16 show an exemplary embodiment of the 
vessel centers processor of FIGURE 5; 

FIGURE 17A schematically shows an exemplary vessel 
structure corresponding to a single vessel center; 
20 FIGURES 17B, 17C, and 17D schematically show the 

recursive eroding of the vessel structure of FIGURE 17A; 

FIGURE 18A schematically shows an exemplary vessel 
structure corresponding to a pair of overlapping vessel 
centers; 

25 FIGURES 18B, IBC, and IBD schematically show the 

recursive eroding of the vessel structure of FIGURE ISA; 

FIGURE 19 schematically shows the decomposition of a 
square moving window for use in recursive erosion into an upper 
left component and a lower right component; 
30 FIGURE 20 shows an exemplary method for the first 

pass of an exemplary two-pass recursive erosion process; 

FIGURE 21 shows an exemplary method for the second 
pass of an exemplary two -pass recursive erosion process; 

FIGURE 22A shows a synthetic structure consisting of 
35 two isolated circular flood-filled regions; 

FIGURE 22B shows the results of the recursive erosion 
method of FIGURES 20 and 21 as applied to the synthetic 
structure of FIGURE 22A; 

- 7 - 
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FIGURE 22C shows a synthetic structure consisting of 
two overlapping circular flood- filled regions; 

FIGURE 22D shows the results of the recursive erosion 
method of FIGURES 20 and 21 as applied to the overlapping 
5 synthetic structure of FIGURE 22C; 

FIGURE 23 shows an exemplary path tracking process 
formed in accordance with an embodiment of the invention that 
employs a representative set of vessel center tags; 

FIGURE 24 shows a robust general tracking system 
10 formed in accordance with an embodiment of the invention that 
employs a representative set of vessel center tags to ensure 
complete tracking of the full vascular skeleton; 

FIGURE 25 shows an exemplary embodiment of the 
direction processor of FIGURE 23; 
15 FIGURE 26A schematically shows the exemplary vessel 

structure of FIGURE 18A with initial geometric contours 
arranged about the vessel centers; 

FIGURE 26B schematically shows the propagating of the 
initial geometric contours of FIGURE 26A; 
20 FIGURE 27 shows an exemplary embodiment of the vessel 

boundaries contour fitting; 

FIGURE 28 shows the exemplary embodiment of the 
contour propagation processor of FIGURE 27; and 

FIGURE 29 shows an exemplary embodiment of a three - 
25 dimensional display processor for displaying the vascular 
skeleton generated by the method of FIGURE 24 , 

Detailed Description of the Preferred Embodimanta 

With reference to FIGURE 4, a magnetic resonance 

30 imaging system that suitably practices angiographic imaging in 
accordance with an embodiment of the invention is described. 
Although the invention is described herein with respect to a 
magnetic resonance imaging embodiment, those skilled in the art 
will appreciate that the invention is applicable to a broad 

35 range of angiographic modalities and techniques, including but 
not limited to contrast -enhanced magnetic resonance 
angiography, non-contrast enhanced magnetic resonance 
angiography, computed tomographic angiography, and fused 
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magnetic resonance/computed tomography angiographic techniques . 
The invention is also suitably practiced in conjunction with 
either white blood angiography (WBA) or black blood angiography 
(BBA) . 

5 With reference to FIGURE 4, a magnetic resonance 

imaging (MRI) scanner 10 typically includes superconducting or 
resistive magnets 12 that create a substantially uniform, 
temporally constant main magnetic field Bq along a z-axis 
through an examination region 14. Although a bore-type magnet 

10 is illustrated in FIGURK 4, the present invention is equally 
applicable to open magnet systems and other known types of MRI 
scanners. The magnets 12 are operated by a main magnetic field 
control 16 . Imaging is conducted by executing a magnetic 
resonance (MR) sequence with the subject being imaged, e.g. a 

15 patient 42 in a magnetic resonance angiography (MRA) session, 
placed at least partially within the examination region 14, 
typically with the region of interest at the isocenter. 

The magnetic resonance sequence entails a series of 
RF and magnetic field gradient pulses that are applied to the 

20 subject to invert or excite magnetic spins, induce magnetic 
resonance, refocus magnetic resonance, manipulate magnetic 
resonance, spatially and otherwise encode the magnetic 
resonance, to saturate spins, and the like. More specifically, 
gradient pulse amplifiers 20 apply current pulses to a whole 

25 body gradient coil assembly 22 to create magnetic field 
gradients along x-, y-, and z-axes of the examination region 
14. 

An RF transmitter 24, preferably digital, applies RP 
pulses or pulse packets to a whole-body RP coil 26 to transmit 

30 RF pulses into the examination region. A typical RF pulse is 
composed of a packet of immediately contiguous pulse segments 
of short duration which taken together with each other and any 
applied gradients achieve a selected magnetic resonance 
manipulation. The RF pulses are used to saturate, excite 

35 resonance, invert magnetization, refocus resonance, or 

- 9 - 
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manipulate resonance in selected portions of the examination 
region . 

For whole-body applications, the resulting resonance 
signals, generated as a result of a selected manipulation, are 
5 also picked up by the whole-body RF coil 26. Alternately, for 
generating RF pulses in limited regions of the subject, local 
RF coils are placed contiguous to the selected region. For 
example, as is known in the art, an insertable head coil 28 is 
inserted surrounding a selected brain region at the isocenter 

10 of the bore. Other surface coils or other such specialized RF 
coils may also be employed. For example, the RF system 
optionally includes a phased array receive coil (not shown) 
which enables partial parallel imaging (PPI) techniques known 
to the art. In one suitable embodiment, the whole-body RF coil 

15 26 induces resonance and the local RF coil or coil array 
receives magnetic resonance signals emanating from the selected 
region. In other embodiments, the local RF coil both excites 
and receives the resulting magnetic resonance signals . 

Regardless of the RF coil configuration and the 

20 application thereof, the resultant RF magnetic resonance 
signals that are picked up by one or another of the RF coils is 
received and demodulated by an RF receiver 32. A sequence 
control processor 34 controls the gradient pulse amplifiers 20, 
the RF transmitter 24, and the RF receiver 32 to produce 

25 integrated MRI pulse sequence and readout waveforms that 
generate the magnetic resonance (MR) signals and optional 
echoes, provide appropriate encoding gradients to spatially 
encode the resultant MR response, and coordinate MR pickup and 
receive operations. 

30 The MRI sequence typically includes a complex series 

of magnetic field gradient pulses and/or sweeps generated by 
the gradient amplifiers 20 which along with selected RF pulses 
generated by RF coils 26, 28 result in magnetic resonance 
echoes that map into k- space. The resultant magnetic resonance 

35 data is stored in a k-space memory 36. The k-space data is 
processed by a reconstruction processor 38, which is typically 

- 10 - 
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an inverse Fourier transform processor or other reconstruction 
processor known to the art, to produce a reconstructed image 
representation that is stored in an image memory 40. £ 
magnetic resonance angiography (MRA) , a patient 42 is imaged by 
5 the MRI system 10 using imaging conditions that provide either 
an enhanced high intensity signal from the vascular regions 
(white blood angiography) or a suppressed low intensity signal 
from the vascular regions (black blood angiography) for the 
vascular regions. The enhanced or suppressed vascular signal 

10 provides contrast for the vascular system in the resultant 
image. In the exemplary embodiment of FIGURE 4, the carotid 
area of the patient 42 is imaged. Optionally, the patient 
receives a magnetic resonance contrast agent 44 to improve the 
vascular contrast, i.e. contrast -enhanced MRA is performed. An 

15 angiographic sequence such as a time of flight (TOF) white 
blood sequence, a gradient recalled black blood sequence which 
can be used in combination with large bipolar gradients or pre- 
saturation RP pulses, a spin-echo (SE) type black blood 
sequence which utilizes wash-out effects and which does not 

20 require large bi-polar gradients or pre-saturation RF pulses, 
or other appropriate angiographic sequence is employed. The 
selected angiographic sequence is applied by the sequence 
control processor 34 and effectuates acquisition of the MRA 
data, 

25 Regardless of the choice of MRA imaging modality, the 

k- space data is collected in the k- space memory 36 and 
reconstructed by the reconstruction processor 38 to generate 
the MRA volume image data 40. The MRA volume image data is in 
the form of a three-dimensional gray scale image representation 

3 0 of the examined area of the patient, which has good contrast 
for the vascular system relative to other body tissues of the 
patient 42. Typically a three dimensional image data 
comprising multiple slices is acquired to provide volume 
imaging of the patient 42. However, it is also contemplated 

35 that only a single image slice is acquired in a patient imaging 
session. 

- 11 - 
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A pre-processor 46 conforms the angiographic data to 
a standardized format preparatory to further processing. For 
example, the pre-processor 46 optionally smooths the data, 
converts the data into an isotropic format, re- sizes the image, 
5 performs pre-f iltering to remove noise or extraneous features 
such as non-vascular contrast that interferes with the vascular 
image, and the like. For BBA, the pre-processor preferably 
intensity- inverts the image so that the vascular areas are 
represented by high intensities. 
10 An edge volume processor 48 receives the output of 

the pre-processor 46, and applies a mathematical transformation 
such as second order spatial differentiation to obtain an edge- 
enhanced volume that particularly emphasizes the edges of the 
vasculature . 

15 A vessels center tagger 50 receives the output of the 

edge volume processor 48 and searches the edge volume for 
vessel centers, e.g. on a slice-by-slice basis. The located 
vessel center points in each plane are tagged. Vessel centers 
corresponding to overlapping vessel images are advantageously 

20 identified and particularly noted. The vessel centers along 
with the overlap tags are supplied along with the edge volume 
image and the pre-processed volume image to a vessel 
segmentation engine 52. 

The segmentation processor 52 processes the 

25 angiographic data to segment the vascular structure. The 
segmenting can be performed in several ways, including: identifying 
the vasculature as a set of voxels corresponding to the imaged 
vascular structures; separating of the vascular regions of the 
image from the non-vascular regions of the image ; assigning a 

30 first value to voxels corresponding to imaged vascular 
structures, and a second value to voxels corresponding to 
imaged non-vaacular structures; and removing contrast 
corresponding to non-vascular structures from the image. 

The segmented vascular information is preferably 

35 graphically displayed on an appropriate user interface 54, typically 
as a two-dimensional or three-dimensional graphical image 

- 12 - 
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representation. The vascular information can of course also be 
stored electronically or on a magnetic storage medium, printed on 
paper, and et cetera {not shown) . 

With reference to FIGURE 5, a method 70 for 
5 post-acquisition processing of acquired raw angiographic volume 
image data Ira„(x,y,z) 72 is described in overview. The raw 
data Ij.a„Cx,y,z) 72 is pre-processed 74 to produce a pre- 
processed angiographic volume I(5(x,y,z) 76 which conforms to a 
pre-selected data format and which has non-vascular contrast 

10 substantially removed therefrom. For example, the pre- 
processor 74 can convert the image format to one formed of 
isotropic voxels, pre-filter the data to remove noise, re-size 
the data, remove extraneous non- vascular contrast, and so 
forth. In the case of black blood angiographic (BBA) data, 

15 the pre-processor 74 preferably removes non-vascular "black" 
contrast such as bone linings, air pockets, and muscle tissues 
that typically are imaged with low intensities similar to the 
black blood. The pre-processor 74 preferably also intensity- 
inverts BBA data so that the vascular regions correspond to 

20 high intensity areas of the pre-processed image Io(x,y,z) 76. 

With continuing reference to FIGURE 5, the 
pre-processed volume Io{x,y,z) 76 is input to a 
three-dimensional high curvature feature processor 78 which 
outputs a three-dimensional edge volume Iedg9Cx,y,z) 80. The 

25 edge volume 80 emphasizes the edges of the vascular structures. 
Processing of angiographic data based on the edge volume 80 is 
advantageous because it retains the vessel lumen information 
and will accurately reflect narrowing of the blood vessels due 
to plaque buildup. 

30 With continuing reference to FIGURE 5 and with 

further reference to FIGURE 6A, a vessel centers processor 82 
receives the edge volume ledge (x, y, z) 80 aiid processes it on a 
slice-by-slice basis to identify vessel centers 84 
substantially throughout the volume of interest which are 

■ - 13 - ■ 

13 



wo 03/041584 



PCT/US02/36048 



representative of the imaged vascular system or systems. 
FIGURE 6A shows exemplary parallel slices 100 arranged 
perpendicular to the paper, with identified vessel centers 
marked or tagged by filled black circles 102. It will be 
5 appreciated that the shown vessel tags 102 are representative 
of a vessel portion having a lower portion 104, a bifurcation 
point 106, and two vessel portions 108, 110 extending upward 
from the bifurcation point 106. Although not immediately 
apparent from the two-dimensional schematic representation of 

10 FIGURE 6A, it will be appreciated that each parallel slice 100 
extends two-dimensionally into and out of the paper and 
typically contains vessel centers across the two-dimensional 
surface, so that the combined slices 100 provide a set of 
vessel center tags 102 which are representative of the complete 

15 three-dimensional vascular system. 

With continuing reference to FIGURES 5 and 6A, the 
vessel center tags 84, 102 preferably include resolved vessel 
centers of overlapping vascular structures. In the exemplary 
schematic FIGURE 6A, the slice 100^ just above the bifurcation 

20 106 includes two vessel centers 102i, 1022 which are likely to 
have overlapping vessel boundaries in the angiographic image. 
The vessel centers processor 82 advantageously recognizes such 
overlapping conditions. 

The vessel center tags 84 form a valuable database of 

25 global information about the vascular structures contained in 
the image. The tags 84 provide a valuable input for many 
contemplated angiographic post -acquisition processing 
applications. In one exemplary application (not shown), the 
centers 84 are plotted to provide a graphical representation of 

30 the vascular system. Given a sufficiently high density of 
vessel centers corresponding to a close spacing of the parallel 
planes 100, a reasonable representation of the vascular paths 
is provided thereby. This paths representation does not, 
however, include the vessel lumen information which is critical 

- 14 - 
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for identifying clinically significant plaque buildup and 
resultant vessel flow constriction which is typically of great 
interest to medical personnel. Neither does this simple 
approach provide for artery-vein separation or other 
5 sophisticated vascular analyses. 

A suitable processing method which incorporates the 
global vessel center tags 84 to effectuate rapid and accurate 
vessel segmentation including accurate vessel lumen 
reconstruction (even in the case of overlapping vessel images) 

10 and optional artery-vein separation is described with 
continuing reference to FIGURES 5 and with further reference to 
FIGURES 6B and 6C. A vessel segmentation processor 86 receives 
the tags 84. The processor 86 selects a starting location 
{x,y,z) 88 in the three-dimensional vascular system. The seed 

15 88 can correspond to one of the tags 84. A direction processor 
90 receives the location 88 and analyzes the three-dimensional 
angiographic volume Io(x,y,z) 76 in the neighborhood of the 
location 88 to deteirmine the vessel direction 92 and the plane 
orthogonal to the vessel direction 94 at the location 88. The 

20 vessel segmentation processor 86 analyzes the orthogonal plane 
94 to determine the vessel boundaries 120 (FIGURE 6B) in the 
plane 94, taking into account tagged vessel overlaps as 
necessary, and also calculates the next vessel center nearest 
to the point 88. The segmentation process is repeated until 

25 the global set of vessel tags 84 are fully visited. With the 
vessel boundaries 120 determined, a three dimensional segmented 
. vasculature 96 can be constructed, e.g. by interpolation 
between the planes 124 and along the vessel centers 126 as 
shown in FIGURE 6C. 

30 By generating and employing the vessel center tags 

84, the robustness of the vascular tracking system 70 is 
greatly improved compared with trackers of the prior art . 
Prior art tracking systems operate in local space, working from 
a starting seed and tracking the direction locally. Such 

- 15 - 
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systems can fail at bifurcations, vessel gaps, vessel overlaps, 
and at other localized image structures which cause the vessel 
direction to become indeterminate. By incorporating the tags 
84, a global database of information about the vasculature is 
5 included into the tracking. The tags 84 overcome the localized 
image structures by providing a global representation of the 
entire vascular system. For example, a vessel gap is 
identifiable by the existence of additional vessel center tags 
84 positioned beyond the gap. Furthermore, the vessel tags 
10 identify overlapping vessels so that the tracker can account 
for the plurality of vessels at a point of overlapping or 
furcation. 

Having provided an overview of the tracking system 
embodiment 70 formed in accordance with the invention with 

15 reference to FIGURES 5 to 6C, the details of the exemplary 
segmentation system 70 are now described. 

With reference to FIGURE 7, a suitable embodiment of 
the pre-processor 74 is described. The raw angiographic volume 
lraw(x,y,z) 72 is input to an isotropic volume processor 140 

20 which converts the gray scale volume image data 72, which may 
or may not be formed of cubic voxels, into an isotropic gray 
scale angiographic image volume 142 which is formed of cubic 
voxels of a selected resolution 144. Of course, if the raw 
data 72 as-acquired is isotropic, the conversion 140 is not 

25 applicable. 

With continuing reference to FIGURE 7, the isotropic 
data is optionally down-sized using an edge -pre serving 
processor 146 which in a suitable embodiment employs a 
shrinking and translating transform such as a discrete wavelet 

30 transform (DWT) 148 for edge-preserving re-sizing, to produce 
a re-sized image 150 which is down-sized by a factor K 152. 
The DWT advantageously substantially preserves the high 
frequency edges of the image so that accurate blood vessel 
lumen information is retained. The wavelet transform is 

35 applied according to: 

- 16 - 
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(1) 



where $ is the isotropic image volume 142, a is a scale factor 
and b is a location/time parameter. The continuous wavelet 
transform is given by convolution as : 



The mother function 3> (herein the isotropic image volume 142) 
should satisfy admissibility criterion to ensure that it is a 
localized zero-mean function. Equation (2) can be discretized 
by restraining a and b to a lattice of the form a=2"^ to get a 
wavelet basis in which the dilations and translations (scaling 
and shifting) of the mother function $ is given by: 



where s and 1 are integers that scale and dilate the mother 
function 4) to generate wavelets, such as those of the 
Daubechies wavelet family. The scale index s indicates the 
wavelet's width and the location index 1 gives its position. 
Those skilled in the art will recognize that by applying 
equation (3) to the isotropic image volume 142, the image is 
rescaled or dilated by K 152 which is a power of 2 , and 
translated by an integer. The image is transformed into four 
down-scaled images, three of which contain horizontal, 
diagonal, and vertical details, respectively, while the fourth 
image contains the desired downsized angiographic volume 150. 



volume, the BBA image is advantageously transformed by a 
pseudo-white-blood-angiographic (pseudo-WBA) processor 156 
which: (i) replaces extraneous black non-vascular contrast with 
an alternative intensity which does not closely correspond with 




(2) . 



(X) = 2 ^^{2-=x-l) 



(3) 



In the case 154 of a black blood angiographic (BBA) 
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the black blood intensity; and (2) intensity- inverts the image 
so that the BBA volume has a "white-blood-type" contrast in 
which the vascular areas appear at high intensity. The 
processor 156 produces a pseudo-WBA image volume 158 which has 
5 the intensity polarity of WBA and is suitable for further 
processing in the same manner as a WBA image volume. 

With continuing reference to FIGURE 7, the image 
volume 150 (for WBA) or the image volume 158 (for BBA) is input 
to an intensity-stretching processor 160 which produces the 
10 pre-processed angiographic volume Io(x,y,z) 76. This processor 
increases the intensity range by increasing the difference 
between the vessels intensity and the background intensity. In 
a suitable embodiment, the intensity-stretching processor 160 
normalizes the voxel intensities according to: 

2W-1 

■^output ~ ^ input ^ ~7 ir~Tr 

■'max -^niin 



15 where N is the number of bits per voxel, is the maximum 

intensity in a slice and is the minimum intensity in a 

slice. The intensity-stretching processor 160 improves the 
sensitivity to the vasculature by raising the vessels intensity 
above the background. It is most effectively used in 

20 conjunction with noise-filtering processors which remove 
isolated, "hanging" noise voxels and otherwise suppress non- 
vascular intensities . 

With reference to FIGURES 8, a suitable embodiment of 
the pseudo-WBA processor 156 is described. The BBA volume 

25 iBBh^^tV,^) ISObbr is processed to remove any non-vascular dark 
or black regions by constructing 170 a binary bone-air-muscle 
mask 172 and assigning 174 a gray intensity level to the non- 
vascular black areas identified by the mask 172. The gray 
intensity level assigned by the assignment processor 174 

30 advantageously corresponds to the average value of the gray 
non-vascular areas which are not masked by the mask 172. Thus, 

-IS- 
IS 
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the black non- vascular regions are replaced by a uniform 
average gray intensity. 

The output of the assignment processor 174 is a 
vessels-tissue gray scale BBA volume lvt(x,y,z) 176 which has 
5 clearly distinguishable vascular regions of very low intensity, 
i.e. black blood, on a higher intensity, gray background. The 
vessels-tissue gray scale BBA volume l^(.x,y,z) 176 
advantageously replaces the non-vascular tissues which 
typically appear black in BBA images by a higher intensity 

10 corresponding to the average intensity of the gray areas of the 
BBA image. Those skilled in the art will recognize that the 
volume Ivt{x,y,z) 176 overcomes the well-known difficulty in BBA 
image interpretation that certain non-vascular tissues, 
particularly air pockets, bone linings, and muscle tissues, 

15 appear black in BBA images and can interfere with, obscure, and 
confuse interpretation of the vascular information contained 
within the BBA image ISObba- 

With continuing reference to FIGURE 8, the vessels- 
tissue gray scale BBA volume lvt(x,y,z) 176 is intensity- 

20 inverted 178 to produce the pseudo-WBA image volume 158 in 
which blood vessels are represented by high intensity, i.e. 
appear white. With the intensity of a voxel of the gray scale 
BBA volume Ivt(x,y,z) 176 designated x^t and the corresponding 
voxel of the pseudo-WBA image volume 158 designated x„ba, the 

25 intensity inversion is according to: 

where n is the number of bits per voxel. For example, if n=8 
then XHBA=255-Xvt, and so a low intensity black pixel having x„t=0 
converts to a high intensity Xkba=255, Although a linear voxel 
intensity inversion transformation has been described with 
30 reference to equation (5) , other intensity inversion 
transformations are also contemplated, such as non-linear 
transformations which advantageously account for non- 
linearities in the intensity data, 
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With reference to FIGURE 9, a suitable embodiment of 
the bone -air-muscle mask processor 170 is described. A slice 
selector 190 selects a slice 192 from the BBA volume ISObba. 
Each pixel in the slice 192 is classified 194. The number of 
5 classes 196, designated herein as K, corresponds to the number 
of distinct types of pixels in the slice 192. As mentioned 
previously, typical BBA images contain two types of pixels: (1) 
black pixels that correspond to vascular structures or to 
certain non-vascular structures such as air pockets, bone 

10 linings, and muscle tissues; and (2) gray pixels that 
correspond to most organs and other tissue types and form much 
of the image background. Thus, the number of classes 196 is 
selected as K=2 , 

With continuing reference to FIGURE 9, the class 

15 assignment processor 194, a suitable embodiment of which will 
be described later, assigns each pixel of the slice 192 as 
either a black pixel or a gray pixel, corresponding to 
bone /air /vascular structures and tissue background 
respectively. The classification results are suitably 

20 expressed as an intermediate binary slice mask 198, The mask 
198 is intermediate because it does not differentiate between 
vascular and non-vascular black pixels, whereas the desired 
mask should identify only the non-vascular black regions. 

It will be appreciated that the vascular regions are 

25 tubular in structure, whereas non-vascular regions are 
typically large annular structures corresponding to bone 
linings or large irregular structures corresponding to air 
pockets. Thus, the mask processor 170 removes the vascular 
regions from the intermediate mask 198 using a mathematical 

30 reduction/ enlargement sequence 200. In a suitable embodiment, 
an erosion process applies a moving window 202 to erode the 
edges of the black regions of the mask by a pre -selected 
amount . The narrow tubular vascular structures are thus 
completely removed whereas the non-vascular structures, which 

35 are typically larger, have only the edges reduced. A 
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subsequent dilation process enlarges the non-vascular areas 
back to their original uneroded size. The resulting slice mask 
204 contains only the non- vascular structures with the blood 
vessels removed, and is suitable for use in identifying the 
.5 non-vascular black regions of the image. 

With continuing reference to FIGURE 9, the slice mask 
generation process is repeated on a slice-by-slice basis 206, 
208 to produce a mask slice 204 corresponding to each slice of 
the BBA volume ISObba- The mask slices 204 are combined 210 to 

10 generate a mask volume 172 . 

With reference to FIGURE 10, a suitable embodiment of 
the class assignment processor 194 is described. The described 
processor 194 embodiment is based on constructing a 
parameterized statistical model including a Gaussian 

15 distribution. The mean and standard deviation of the Gaussian 
distribution are the optimized parameters. The class 
assignment processor -194 is described with reference to N 
pixels in the image slice 192, and with respect to the K 
classes 196. Although K=2 for typical two-class BBA data, 

20 inclusion of additional classes is also contemplated for 
specific imaging situations, and so the class assignment 
processor 194 is described for an arbitrary number of classes 
K 196. The pixels are designated herein by i where lii^N, and 
the classes are designated by k where l^k^K. The pixel 

25 intensities are expressed herein as x^. The pixel 
classifications are expressed as yi,,;, which indicates the 
statistical probability that the pixel i is a member of the 
class k. A Gaussian probability model 220 is used to calculate 
unweighted pixel probabilities 222 according to: 



30 where: Pi,jc=P (Xi|yi, <E>) is the unweighted Bayesian probability 
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that a pixel of intensity Xi is in the class k; and Ou are 
the mean and standard deviation of the Gaussian probability 



and D is a dimensional parameter, here D=l. 

With continuing reference to FIGURE 10, the described 
embodiment of the class assignment processor 194 optimizes the 
classification probabilities by maximizing the a posteriori 
probability by iteratively optimizing the Gaussian model 
parameters ^=[1x^,0^,] and a weighting w,,. The pixel 

classification values yi^^ and the distribution parameters 
4'={/ik/ crj are assigned initial values 224, and an initial 
probability weighting Wj, is computed. Suitable initial values 
are given by: 



where ? is the covariance and the superscript "o" indicates the 
zeroth iteration. The initial mean and standard deviation 224 
are copied 226 into the past-iteration or old mean and standard 
deviation memoiry 228. Using the Gaussian probability model 
220, the unweighted pixel probability Pi,^ 222 is calculated 
using the old mean and standard deviation 228, and an updating 
of the weighted classification values yi,^"'^^ is performed 230 
according to: 



•distribution model; ^ are the model parameters, e.g. ^={^^.,o^;}; 




(7) 





(8), 



where the superscript n indicates the current iteration, the 
superscript n+1 indicates the next iteration, and Wi^" is a 
weighting given by: 
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(9) . 



The updated classification probabilities yi.ic"*^ 232, which are 
weighted probabilities, are thus obtained. An iterative 
updating of the statistical model parameters is also performed 
234 according to: 



Eyi%(^i-Pt)^ 
£yi% 



to generate updated statistical model parameters 236, 

With continuing reference to FIGURE 10, the iterative 
process is terminated using a suitable convergence criteria. 
For example, convergence can be indicated by computing 238 a 
10 difference between the old mean fi^^" and the new mean /ifc"*^ for 
each class k and checking 240 for convergence by comparing the 
differences with pre-set criteria. The iterative process can 
also optionally be terminated after a pre-set number of 
iterations. 

15 Once the iterative process is complete, an assignment 

optimizer 242 assigns the optimized classification to each 
pixel of the intermediate slice mask 198. For each pixel, the 
final weighted probabilities yi,,, of each class for a pixel are 
obtained and the class k that has the highest corresponding 

20 weighted probability yi,K for that pixel is selected as the 
classification of that pixel. 

A suitable embodiment of an exemplary pre-processor 
74 has been described with reference to FIGURES 7-10. The 
exemplary pre-processor 74 generates an isotropic, edge- 

25 preserving, optionally re-sized and intensity-inverted image 
Ic>{x,y, z) 76 which has WBA intensity polarity. Of course, 
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additional or other pre-processing steps are also contemplated, 
auch as additional noise filtering, smoothing, etc. 

With reference returning to FIGURE 5, the pre- 
processed angiographic volume Io(x,y,z) 76 is input to the 
5 three-dimensional high curvature feature processor 78 which 
outputs an edge volume ledge (x,y, z) 80. 

With reference to FIGURE 11, a suitable embodiment of 
the edge processor 78 is described. The edge processor 78 
approximates differentiation of the angiographic volume 

10 Io(x,y,z) 76 by convolving the image 76 with a second order 
Gaussian derivative. A voxel 262 is chosen 260. A Gaussian 
design processor 264 calculates a Gaussian and its derivatives 
2 66 based on initial coefficients 268. The calculation of the 
Gaussian derivatives preferably takes advantage of known 

15 closed-form expressions for the Gaussian derivatives. Although 
the convolving is ordinarily a calculation of order O(k^) where 
k is the number of number of voxels forming the Gaussian kernel 
266, by convolving in separable mode the calculation is reduced 
to order 0(3k) . Thus, the convolving reduces to: a 

20 convolution in the x-direction 270 which produces an x- 
intermediate result 272; a convolution in the y-direction 274 
which produces a y- intermediate result 276; and a convolution' 
in the z-direction 278 which produces the final result 2 80. 
The gradient processor 282 combines the one -dimensional 

25 derivatives estimated by the convolving to form an edge volume 
voxel according to: 




(11) 



where I is the pre-processed angiographic volume Io{x,y,z) 76 
and the partial derivatives are evaluated at the location 
(x,y,z) of the selected voxel 262. After repeating 284 for 
30 every voxel, the edge volume I,dge(x,y, z) 80 is generated. 
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It will be appreciated that the Gaussian is a 
function with a variance o. If a is large, the Gaussian 
convolution can take prohibitively long, especially for higher 
order derivatives. For calculation of the edge volume 80, 
5 second order derivatives are used. In a preferred embodiment 
which greatly improves calculation speed, the convolving is 
performed using a z-transform decomposition of the Gaussian 
into sinusoidal and cosinusoidal components. 

With reference to FIGURE 12, a suitable embodiment 

10 for performing the z-transform decomposition is described. For 
a scale a 2 90 and a one-dimensional coordinate x 292, a cosine 
processor 294 calculates a cosine term 296, and similarly a 
sine processor 298 calculates a sine term 300. The cosine and 
sine terms are multiplied 302, 304 by the cosine and sine z- 

15 transform coefficients 306, 308, respectively, to produce the 
cosine and sine z-transform components 310, 312, respectively. 
These components 310, 312 are additively combined 314 to 
produce the z-transform 316. 

With reference to FIGURE 13, an exemplary time-of- 

20 flight white blood angiographic image 330 of the carotid area 
is shown. The image 330 includes five isolated vessels 332 as 
well as one pair of overlapping vessels 334. 

With reference to FIGURE 14, the edge volume 340 
corresponding to the image 330 is shown. The edge volume is 

25 formed according to the method described with reference to 
FIGURES 11 and 12. Each of the isolated vessels 332 is shown 
to have a well defined closed contour 342 corresponding 
thereto. Similarly, the overlapping vessels pair 334 has a 
corresponding well defined closed contour 344. Besides the 

30 contours 342, 344, the remainder of the edge volume 340 is 
substantially free of features, except for a few low intensity 
enclosed regions 346 associated with non-vascular features. 
These non-vascular edges 346 are lower intensity versus the 

- 25 - 

25 , 



wo 03/041584 



PCTrtJS02/36048 



vessel contours 342, 344. 

With reference back to FIGURE 5, the vessel centers 
processor 82 receives the edge volume 80 and identifies a 
plurality of vessel center tags 84 representative of the 
5 vascular system. The vessel center tags 84 identify the vessel 
centers in three-dimensional space, e.g. by coordinates 
(x,y,z), and also advantageously identify vessel centers 
corresponding to overlapping vessel images. 

With reference to FIGURE 15, a suitable embodiment of 

10 the vessel centers processor 82 is described. The approach of 
this embodiment is to divide 370 the three-dimensional edge 
volume 80 into two-dimensional edge volume slices 372. Each 
slice 376 is selected 374 for processing in turn. The edge 
volume slice 376 is converted 378 into a binary slice 3 80 by 

15 setting the non-zero values to one 382. The binarization 378 
can be effectuated, for example, by using a very low 
thresholding value. As seen in the exemplary edge volume slice 
340 of FIGURE 14, the edge volume image elements are 
substantially zero (i.e., black} except at the edge contours 

20 342, 344, and so a low thresholding value will appropriately 
set the vessel contours 342, 344 to binary one imposed on a 
binary zero background, thus forming closed binary contour 
edges of the vessel structures . The vessel structures are 
flood-filled 384 to form a slice having flood-filled vessel 

25 structures 386. There are a number of well-known methods for 
flood- filling closed contours, and one or more appropriate 
algorithms are typically available as standard functions in 
signal processing software development environments. The flood- 
filled regions are counted 388 to produce a count and location 

30 390 of the vessel structures in the slice 386, In a suitable 
embodiment, a run length encoder method performs the counting. 
Run length encoders are well-known algorithms which are 
commonly employed in JPEG image compression, and functions 
which implement a run length encoder or a variation thereof are 

35 often available in signal processing software toolboxes or 
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libraries. 

With continuing reference to FIGURE 15 and with 
further reference to FIGURE 16, a vessel centers tagger 400 
receives the slice 386 with the flood filled vessel structures 
5 and the count and location of the vessel structures 390. For 
each flood-filled vessel structure of the slice 386, the vessel 
centers tagger 400 identifies one or more vessel centers 402 
corresponding thereto. The vessel centers tagger 400 is 
suitably embodied by a recursive erosion process to be 

10 described later. Advantageously, the vessel centers processor 
82 also particularly tags overlapping vessel structures by 
identifying 404 vessel structures having more than one vessel 
center 402, since a plurality of vessel centers 402 
corresponding to a single vessel structure indicates the 

15 singular vessel structure is formed of overlapping vessel 
images. Appropriate vessel overlap tags 406 are generated by 
the vessel overlap identifier 404. By cycling 408 through the 
slices 372 a set of vessel centers and overlap tags 84 is 
generated that is representative of the imaged vascular system. 

20 With reference to FIGURES 17A to 17D, a suitable 

method for implementing the vessel centers tagger 400 is 
described. A flood filled vessel structure 420 is successively 
eroded as shown in FIGURES 17B through 17D using a moving 
window 422. The recursive erosion process iteratively removes 

25 outer portions 424 of the vessel structure 420, leaving 
successively smaller vessel structure portions 426, 428. The 
final remaining vessel structure portion 430 after the 
recursive eroding tagged as the vessel center, e.g. its (x,y,z) 
spatial coordinates are recorded in the table of vessel centers 

30 84. It will be recognized that the flood filled vessel 
structure 420 of FIGURE 17A erodes to a single vessel center 
430. 

With reference to FIGURES ISA to 18D, the method of 
FIGURES 17A to 17D is applied to a flood filled vessel 
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Structure 440 which images overlapping vessels. During the 
initial stages of the recursive eroding, the moving window 422 
removes outer layers 444 leaving a single reduced vessel 
structure portion 446. However, as the iterative eroding 
5 continues, two distinct vessel structure portions 448, 450 
become separately distinguishable. Ultimately, the vessel 
structure portions 448, 450 erode down to identifiable vessel 
centers 452, 454. Since the vessel centers 452, 454 correspond 
to a single flood filled vessel structure 440, the vessel 

10 centers processor 82 identifies 404 the flood-filled image 440 
formed of overlapping vessel images corresponding to the vessel 
centers 452, 454. 

The moving window 422 shown in exemplary FIGURES 17B, 
17C,' 18B, and 18C is a filled circular element. With reference 

15 to FIGURE 19, an exemplary square window 470 is shown which is 
formed from nine pixel elements arranged about a center (r,c) . 
In order to in^rove the speed of the recursive erosion, the 
erosion is advantageously separated into two passes. A first 
pass runs across the slice from a first corner, e.g. the upper 

20 left corner, to a second opposite corner, e.g. the lower right 
cornfer, and uses only the upper left portion 472 of the moving 
window 470. The first pass is followed by a second pass that 
runs across the slice from the second corner, e.g. the lower 
right corner, to the first corner, e.g. the upper left comer, 

25 and uses only the lower right portion 474 of the moving window 
470. 

With reference to FIGURE 20, the first pass of the 
decomposed recursive erosion embodiment of the vessel centers 
tagger 400 is described. For each pixel (r,c} of the slice 386 
30 having flood- filled vessel structures, the top four values A, 
B, C, D around (r,c) are obtained 500, namely: A=Cr-l, c-l) ; 
B=(r-l, c) ; C=(r-1, c+1) ; and D=(r, c-l). It will be 
appreciated that these top four values 502 correspond to the 
upper left portion 472 of the moving window 470 of FIGURE 19. 

- 28 - 
28 



wo 03/041584 



PCT/US02/36048 



A minimum processor 504 selects the minimum value 506 of these 
four values. And add/replacement processor 508 adds one to the 
minimum value 506 and replaces the pixel (r,c) with the new 
value 510. Thus, the transform for the first pass is: 

(r,c) = min {A, Br CD) + 1 (12) 



5 where A, B, C, D are the pixels defined above. The transform 
of equation (12) is repeated 512, 514 for each pixel to produce 
an intermediate LRTB image 516 . 

With reference to FIGURE 21, the second pass of the 
decomposed recursive erosion embodiment of the vessel centers 

10 tagger 400 is described. For each pixel tr,c) of the 
intermediate LRTB image 516, the bottom five values around 
(r,c) are obtained 530, namely: E=(r, c) ; F={r+1, c+1) ; G= (r, 
c+1) ; H=(r+1, c) ; and I=(r+1, c-1) . It will be appreciated 
that these bottom five values 532 correspond to the lower right 

15 portion 474 of the moving window 470 of FIGURE 19. An 
add/ replacement processor 534 adds one to each of the five 
values E, F, G, H, I to create five new values 536. A minimum 
processor 538 selects the minimum value as the new value 540 at 
the pixel (r,c} . Thus, the transform for the second pass is: 
(r,c)^^^ = min{£+l,F+l,G+l,Jf+l, 1+1} (13) 



20 where E, F, G, H, I are the pixels defined above. The 
transform of equation (13) is repeated 542, 544 for each pixel 
to produce the eroded image 546 resulting from an erosion pass. 
The recursive erosion process of FIGURES 20 and 21 is repeated 
until the final vessel centers are obtained. 

25 With reference to FIGURES 22A through 22D, exemplary 

results for recursive erosion according to the process of 
FIGURES 20 and 21 are presented. In FIGURE 21A, a synthetic 
structure 570 is formed of two fully distinct flood-filled 
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circular contours 572, 574. As shown in FIGURE 22B, the 
recursive erosion process performed according to FIGURES 20 and 
21 produces two eroded structures 576, 578 that have clearly 
defined high- intensity centers corresponding to the centers of 
5 the two flood-filled circular contours 572, 574. 

With reference to FIGURES 22C through 2 2D, a 
synthetic structure 580 is formed of two fully distinct flood- 
filled circular contours 582, 584. As shown in FIGURE 2 2D, the 
recursive erosion process performed according to FIGURES 20 and 

10 21 produces two eroded structures 586, 588 that have clearly 
defined high- intensity centers corresponding to the centers of 
the two flood-filled circular contours 582, 584. A high 
intensity connecting line 590 connects the high intensity 
centers 586, 588 but does not detract from the detectability of 

15 the centers 586, 588 because they are positioned at the 
. endpoints of the line 590. 

With reference to FIGURE 23, the exemplary vessel 
segmentation processor described in brief overview with 
reference to FIGURE 5 is described in greater detail. The 

20 segmentation process of FIGURE 23 uses a path tracking process 
600. A starting point 610 lying within the vascular system is 
selected. The starting point is typically a root of the venous 
or arterial system of interest, and can optionally be selected 
manually or identified using an automated system, e.g. by 

25 selecting the vessel structure having the largest area from the 
count and location of vessel structures 390 (FIGURE 15) . The 
starting point is optionally selected from the table of vessel 
center tags 84. The tracking point at (x,y, z) 88 is selected 
612, initially corresponding to the starting point 610. 

30 Tracking from the point (x,y,z) 88 includes (1) 

finding the vessel direction; (2) finding the vessel center; 
(3) finding the extent of the vessel about the vessel center in 
the plane orthogonal to the vessel direction; (4) locating the 
next vessel center; and (5} repeating the process. The 
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direction processor 90 finds the vessel direction 92 and the 
orthogonal plane 94. A vessel centers search processor 612 
searches the vessel centers tags 84 to locate the nearest 
vessel center 614 in space to the point (x,y,z) 88. (Of 
5 course, if the point 88 is selected from the vessel centers 
tags 84, the searching 612 can be omitted) . The segmentation 
processor 86 analyzes the orthogonal plane 94 to determine the 
extent of the vessel about the vessel center 614. The 
boundaries of the vessel along with the vessel center 612 form 

10 a vessel skeleton portion 616. The process is repeated 618, by 
selecting 612 a new point 614. The selecting 612 
advantageously moves a pre -selected distance along the vessel 
direction 92, e.g. a distance equivalent to one voxel for the 
best, tracking resolution. Once the vessel termination is 

15 reached, a vascular path 96 « is obtained. 

The tracking system described with reference to 
FIGURE 23 applies to a vascular path 96'. At bifurcation 
points, the overlap tags 406 (FIGURE IS) are advantageously 
used to identify new starting points 610 at which the path 

20 tracking 600 of FIGURE 23 is repeated. 

With reference to FIGURE 24, a robust general 
tracking system 630 that uses the tags 84 to ensure tracking of 
all the vascular branches is described. The vessel starting 
point 610 is selected 632 from the vessel center tags 84. The 

25 tracking system 600 of FIGURE 23 is applied to generate the 
tracked path 96' corresponding to the starting point 610. 
During the tracking 600, each tracked vessel center 614 is 
marked as visited. After the tracking 610, a check 634 is 
performed to determine whether all the vessel tags 84 have been 

30 marked as visited. If there are unvisited tags, one of the 
unvisited tags is selected 632 as the new starting point 610, 
and the tracking 600 is repeated. Once all the vessel centers 
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84 have been visited, the tracked paths 96" are combined to 
form the full vascular skeleton 96. Because the vessel center 
tags 84 are representative of the imaged vascular system, the 
system is robust. Bifurcation points, tortuous or occluded 
vessels, vessel overlaps, intertwined vessels, partial volume 
averaging and other imaging artifacts, and vessel gaps which 
act alone or in various combinations to produce localized 
disjoints or multiplicities of the vascular path do not prevent 
the general tracking system 630 from tracking all vessel 
branches, because visitation of every vessel branch is ensured 
by visiting every one of the vessel tags 84 . 

With reference to FIGURE 25, a suitable embodiment of 
the direction processor 90 is described. The concepts of 
differential geometry, curves and surfaces are preferably used 
for estimating the direction 92 and the orthogonal plane 94. 
In 2-D images, lines or curves have two directions, one 
direction is along the line or curve and second direction that 
is orthogonal to the line or curve. These directions are the 
vectors corresponding to the eigenvalues of the Hessian matrix 
at a given point. For a 2-D image, the two eigenvalues are the 
two curvatures at that point, i.e. the maximum curvature (kj 
and the minimum curvature (ka) . Using the differential 
geometry, one can represent the amplitude of the curvature and 
its corresponding direction by the eigenvalues and eigenvectors 
of the Weingarten matrix W = Fi'^F;, where Fi is a function of 
the partial derivatives of the image in x and y directions, and 
F2 is the function of the second partial derivatives in x and 
y directions. 

The three dimensional pre-processed angiographic 
volume Io(x,y,z) 76 is modeled as a four- dimensional surface 
given by: 



I(x,y,z) 
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where the first three elements are the spatial coordinates 
(x,y,z) and the fourth element is the image intensity I(x,y,z) 
at the voxel location (x,y,2) . In this case, there are three 
principal curvatures. Two curvatures correspond to the two 
5 orthogonal directions in the cross-sectional plane of the 
tubular blood vessel, while the third principal curvature 
corresponds to the vessel direction or vessel orientation. The 
three directions can be computed using the Weingarten matrix 
which is a combination of the first and second form of the 

10 hypersurf ace, i.e. W = Pr'"F2, where Fi is the fundamental form 
of hypersurf ace and a function of I,, ly and 1^, the three 
partial derivatives of the image volume, and is the second 
fundamental form of hypersurface and is a combination of the 
second partial derivatives I^^, I^y, I^^, lyy, lyj, and 1^^. More 

15 explicitly, 

W = Fi"^ 
wiiere 

F, = I^I^ l^ll I^I^ 

^Jz ^y^z l*^z. 



A Weingarten (W) matrix 652 is constructed 650 according to 
equation (15) . The eigenvectors and eigenvalues 656 of this 
matrix are obtained by diagonalizing 654 the Weingarten matrix, 
such eigenvalue decomposition being well known to those of 

20 ordinary skill in the art. The vessel direction 92 and the 
orthogonal plane 94 are identified 658 from the eigenvalues and 
eigenvectors 656 as follows. The eigenvector direction 
corresponding to the minimum curvature is the vessel direction 
92. The plane passing though the point (x,y, z) whose normal is 

25 the vessel direction 92 is the orthogonal plane 94. 

The orthogonal plane 94 is in general an oblique 
plane relative to the principle axes of the isotropic MRA 
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volume 76. That is, the plane 94 is not in general parallel to 
one of the axial, sagittal, and coronal planes along which MRA 
data is typically measured. Trilinear interpolation of the MRA 
volume 76 is preferably performed to obtain the orthogonal 
5 plane 94 with isotropic pixels. 

With reference returning to FIGURE 23, the vessel 
segmentation processor 86 receives the orthogonal plane 94 and 
estimates the vessel boundaries in the orthogonal plane 94 
which form part of the vessel skeleton 616. 

10 With reference to FIGURE 26A and 26B, a suitable 

embodiment of a method for estimating the vessel boundaries is 
schematically described. FIGURES 26A and 26B schematically 
show estimation of the vessel boundaries for the flood- filled 
overlapping vessel structure 440 of FIGURE 18A. It is known 

15 from the table of vessel centers and overlap tags 84 that there 
are two vessel centers 452, 454 associated with the overlapping 
vessel structure 440 {see FIGURE 17D) . Initial geometric 
contours 680, 682 are arranged around the vessel centers as 
shown in FIGURE 26A. The contours 680, 682 are level set 

2 0 contours defined by the partial differential equation: 

= (SK + V^)\V^>\ - V^^,V(p (16) 

where is the level set contour 680, 682, ek, Vp, and Ve^t are 
speed functions known as the curvature, regional, and gradient 
speed functions, respectively, that control the propagation or 
flow of the contour. The level set contours 680, 682 are 

25 iteratively expanded from their initial position based on the 
speed functions according to equation (16) , and as shown in 
FIGURE 26B, to produces expanded level set contours 684, 686 
which iteratively fill and define the vessel boundaries. The 
region-based level set approach has been previously used in 

30 brain tissue segmentation applications. Its application herein 
to vascular boundary estimation provides improved accuracy and 
speed. Methods previously used for vessel boundary estimation, 
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particularly parametric contour methods, are inaccurate at 
sharp corners and have a tendency to bleed through gaps in the 
vessel structure. 

With reference to FIGURE 27, the level set approach 
5 as applied to vascular boundary reconstruction is described. 
An initial contour 702 is constructed 700 about the vessel 
center selected from the vessel center tags 84 in the 
orthogonal plane 94. The vessel contour is propagated, i.e. 
expanded 704. The propagating 704 is iterative, occurs within 

10 the plane 94, and is constrained by appropriate constraints 706 
such as the distance from the vessel center and contact between 
the propagating 704 contour and neighboring contours. Another 
constraint is the outer boundary 708 of the vessel structure in 
the image. These constraints serve as stopping criteria which 

15 prevent the contour from expanding too far, e.g. beyond the 
vessel structure or into another vessel. The final propagated 
contour defines the vessel boundary 710. The level set 
approach is repeated 712 for each vessel center of the vascular 
structure. For the exemplary case of the vessel structure 440 

20 (FIGURE 18A) the level set approach is repeated 712 for the two 
vessel centers 452, 454 as shown in FIGURES 26A and 26E. 

With reference to FIGURE 28, a suitable embodiment of 
the contour propagator 704 is described. The propagator 704 of 
FIGURE 28 is described with reference to a pair of overlapping 

25 vessels such as are shown in FIGURES 26A and 26B. Those 
skilled in the art can easily make the necessary changes to 
accommodate a vessel structure having only a single vessel 
center associated therewith. Thus, there are two initial 
contours 702, for which initial level set distance map fields 

30 732 are computed 730 based on: the arrangement of the initial 
contours 702 about the vessel centers; the orientation of the 
orthogonal image plane 94; and the selected width of a narrow 
band 734 within which the level set field 732 is applied. 

With continuing reference to FIGURE 28, a new field 
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738 is calculated 736 by solving equation (16) iteratively to 
yield: 

<p"'i = <y - At {V^^(x,y) + V^„d{x,y) -V^„^(x,y) } (17) 

where an integrated speed input 740 includes the three velocity 
terms Vreg(x,y) , Vgrad{x,y), andV5„(x,y). Typically, At=l, i.e. 
each iteration corresponds to a discrete time interval. 
Vreg(x,y) is a regional speed or velocity given by: 

V^^g{x,y) = max{l^p, 0} min{Vp, 0} V- 

where 

a. (18), 

VJx,y) - ■ 



Y[l-2u{x,y)] 

Vgrad{x,y) is a gradient speed or velocity given by: 
V,,(x,y) =V,,,(x,y) -V,,,,,,{K,y) 
where 

max{p"i,x,y) ,0] ^"''(x, y) +niin{ g"{x, y) , 0} D**(x,y) , 
max{g"(x,y),0} D-^{x, y) +inin(p " (x,y) , 0} D*y{x,y) 



(19), 



and Vgraci(x»y) is a curvature speed or velocity given by: 

V^„^(x,y) = 2K"(x,y) [(D°''(x,y))2 + (D°^{x,y))2]i/2 (20). 

In the above equations, y is a damping coefficient, u(x,y) is 
10 a fuzzy pixel membership function whose value lies in the range 
[0,1], Mr is a regional weighting, p" and q" are the x- and y- 
cotnponents of the gradient strength and are calculated from the 
edge volume 80, K"(x,y) is the curvature at location (x,y) , and 
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Vx*, Vy*, Vjt", Vy" are the forward and backward level set 
gradients in the x and y directions which can be expressed in 
terms of finite difference operators of the form D*^'*"'*'* (x,y) 
for forward and backward differences, and as d"'"'''^* (x, y) for 
5 averaged differences. The regional velocity V^eg relates to the 
fuzzy meitibership function u(x,y) which for example could be 
computed using the fuzzy C mean algorithm known to the art . 
The gradient velocity Vg„d relates to the gray scale intensity 
gradients, and the curvature velocity relates to the 

10 curvature in level set space. 

With continuing reference to FIGURE 28, once the new 
level set field calculated according to equation (17) , the 
contour is updated and the level set field is reinitialized 742 
using isocontouring 746 in which the new contour is defined 

15 essentially corresponding to the zero contour of the level set 
field 738. The updating can employ the fast marching method or 
other techniques known to the art. The calculating 742 of the 
updated contour and field 744 is also constrained by selected 
constraints 706 as discussed previously with reference to 

20 FIGURE 27. Of course, since FIGURE 28 refers to an exemplary 
propagation of a pair of overlapping vessels such as are shown 
in FIGURES 26A and 26B, there are two contours 684, 686 (FIGURE 
26B) being calculated, one each about the vessel centers 452, 
454 (FIGURE 18D) . Thus, an appropriate constraint in this 

25 exemplary case is that the two contours 684, 686 not overlap. 
Similarly, if there were three or more contours corresponding 
to an overlapping vessel structure, an appropriate constraint- 
would be the non-overlapping of the three-or more contours. 

With continuing reference to FIGURE 28, the iterative 

30 level set propagating is repeated 748, 750, 752 a plurality of 
times until the contours converge to the two zero-level 
contours 710 corresponding to the two overlapping vessels. 

Although a contour propagation employing a geometric 
contour in region -based level set space including a regional 

3 5 fuzzy membership velocity component is described herein, other 
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methods for finding the vessel boundaries are also 
contemplated, such as methods using parametric contours, 
geometric contours employing other types of distance maps, and 
the like. 

5 With reference returning to FIGURE 24, the exemplary 

robust vessel tracking system 630 generates a full vascular 
skeleton segmentation 96. This extracted information can be 
conveyed to medical personnel in a variety of ways, including 
three-dimensional graphical displays, two-dimensional 

10 projections, selected vascular sub-system displays (e.g., 
limited to the arteriogram) , and the like. 

With continuing reference to FIGURE 24 and with 
further reference to FIGURE 29, an exemplary embodiment of a 
method for generating a suitable three-dimensional display is 

15 described. The vascular skeleton 96 is supplied as an input. 
The tracking system 630 generates a skeleton 96 formed of a 
sequence of vessel contours or boundaries 710 (FIGURES 27 and 
28) lying in two dimensional slices corresponding to the 
orthogonal planes 94 (FIGURE 23) . Thus, the first step is to 

20 operatively interpolate across the boundaries 710. In the 
exemplary embodiment of FIGURE 29, shaded wire mesh method of 
a type known to the art ia employed. A triangulation method 
770 generates a three-dimensional wire mesh of triangular 
elements 772. The normals of the vertex triangles 776 which 

25 determine the solid vascular shape are computed 774. A display 
processor 778 converts the wire mesh into a shaded three- 
dimensional display 780. Optionally, selected graphical 
animation is produced 778 as well, such aa graphical 
highlighting of a selected vascular system, e.g. the 

30 arteriogram. Of course, the three-dimensional display just 
described is exemplary only, and those skilled in the art will 
recognize that many other methods of conveying the segmented 
vascular tree information can also be used. 
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Having thus described the preferred embodiments, the 
invention is now claimed to be: 

1. A method for characterizing a vascular system in a 
three-dimensional angiographic image comprised of voxels, the 
method comprising: 

extracting (374) from the angiographic image a 
two-dimensional slice (376) formed of pixels; 

. identifying (384, 500, 504, 508, 514, 530, 534, 538, 544) 
vessel center points (402) in the two-dimensional slice (376) ; 

repeating (408) the extracting (374) and identifying (384, 
500, 504, 508, 514, 530, 534, 538, 544) for a plurality of 
slices to generate a plurality of vessel center points (84) 
that are representative of the vascular system; and 

segmenting, tracking, extracting, enhancing, or 
identifying (630) vascular information (96) contained in the 
angiographic image using the vessel center points (84) as 
operative inputs. 

2. The method as set forth in claim 1, wherein the 
identifying (384, 500, 504, 508, 514, 530, 534, 538, 544) of 
vessel center points (402) includes: 

flood-filling (384) imaged vascular structures in the 
slice (376) to form filled regions defined by pixels having a 
first value; and 

iteratively eroding (500, 504, 508, 514, 530, 534, 538, 
544) the edges of the filled regions to identify vessel center 
points (402) . 

3. The method as set forth in claim 2 further including: 
prior to the flood-filling (384) , enhancing (78) vessel 

edges by second order spatial differentiation of the 
angiographic image (76) . 
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4. The method as set forth in claim 2, further 
including: 

prior to the flood- filling (384) , enhancing (78) the 
vessel intensity contours by convolving (270, 272, 274) the 
angiographic image (76) with a kernel (266) formed from a 
second or higher order derivative of a Gaussian function. 

5. The method as set forth in claim 4 wherein the 
convolving (270, 272, 274) of the angiographic image (76) with 
a kernel (266) includes: 

decomposing (294, 298, 302, 304, 314) the kernel into 
sinusoidal components (316) ; and 

convolving (270, 272, 274) the angiographic image with the 
sinusoidal components of the kernel. 

6. The method as set forth in any one of claims 2-5, 
wherein the iterative eroding (500, 504, 508, 514, 530, 534, 
538, 544) of the edges of the filled regions includes: 

performing a first erosion pass (500, 504, 508, 514) in a 
first direction using a first moving window (472) ; 

performing a second erosion pass (530, 534, 538, 544) in 
a second direction using a second moving window (474) ; and 

repeating the first and second erosion passes (500, 504, 
508, 514, 530, 534, 538, 544) a plurality of times such that 
the remaining at least one region (452, 454) corresponds to 
vessel center points. 

7. The method as set forth in any one of claims 2-6, 
further including: 

tagging (404) vessel overlaps and vessel furcations (406) 
identified as a plurality of vessel centers (402) corresponding 
to a single filled region. 
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8. The method as set forth in any one of claims 1-7 
wherein the segmenting, tracking, extracting, enhancing, or 
identifying (630) of vascular information (96) includes: 

selecting (612) a first vessel center point (610) ; 

finding (90) a vessel direction (92) corresponding to the 
first vessel center point (610) based on analysis of the 
angiographic image in the three-dimensional neighborhood of the 
first vessel center point (610) ; 

defining a plane (94) in the angiographic image which is 
perpendicular to the vessel direction (92) and contains the 
first vessel center point (610) ; 

estimating (86) vessel boundaries (710) corresponding to 
the first vessel center point (610) in the defined plane (94) ; 

repeating (618, 634) the selecting (612), finding (90), 
defining a plane (94) , and estimating (86) for the plurality of 
vessel center points (84) ; and 

interpolating (770, 774) the estimated vessel boundaries 
(710) to produce a vascular representation. 

9. The method as set forth in claim 8 wherein the 
estimating (86) of vessel boundaries (710) includes: 

defining (700) an initial geometric contour (702) arranged 
about the vessel center (610) and lying in the defined plane 
(94); and 

iteratively optimizing (704) the geometric contour 
constrained to lie in the defined plane (94) . 

10. The method as set forth in claim 9 wherein the 
iterative optimizing (704) of the geometric contour uses a 
level set framework (730, 736, 742) . 



11. The method as set forth in either one of claim 9 and 
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claim 10, wherein the iterative optimizing (704) includes: 

computing (742) a new contour (744) based on a current 
contour (732) and a fuzzy memberahip classification of the 
pixels in the neighborhood of the current contour (732) . 

12. The method as set forth in any one of claims 9-11, 
wherein the iterative optimizing (704) of the geometric contour 
further includes constraining (706) the iterative optimizing 
(704) by at least one of: 

edges (708) of a vascular structure image containing the 

first vessel center (610) ; 

a neighboring vessel boundary (684, 686) ; and 

a pre -determined distance from the vessel center (610) 

about which the geometric contour is arranged. 

13. The method as set forth in any one of claims 8-12, 
wherein the finding (90) of a vessel direction (92) includes: 

constructing (650) a Weingarten matrix (652) ; 

obtaining (654) a plurality of directions (656) by 
eigenvalue decomposition of the Weingarten matrix (652) and 

selecting (658) the vessel direction (92) from the 
plurality of directions (656) , 

14. The method as set forth in any one of claims 1-13, 
further including: 

conditional (154) upon the angiographic image (150) being 
a black blood angiographic image, inverting (156, 160) the 
intensities of the image elements to generate an intensity- 
inverted image (76) . . 

15. An apparatus for characterizing a vascular system, 
the apparatus including an imaging means (10) for acquiring 
imaging data from at least a portion of a patient (42), and a 
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reconstruction means (38) for reconstructing the acquired 
imaging data to generate a three-dimensional image 
representation (72) comprised of voxels and including enhanced 
vascular contrast, the apparatus further comprising: 

a means for extracting (374) from the image representation 
a two-dimensional slice formed of pixels; 

a means for identifying (384, 500, 504, 508, 514, 530, 
534, 538, 544) vessel center points (402) in the two- 
dimensional slice (376) ; 

a means for iteratively invoking (408) the extracting 
(374) and identifying (384, 500, 504, 508, 514, 530, 534, 538, 
544) means for a plurality of slices to generate a plurality of 
vessel center points (84) that are representative of the 
vascular system; and 

a means for segmenting, tracking, extracting, enhancing, 
or identifying (630) vascular information (96) contained in the 
angiographic image using the vessel center points (84) as 
operative inputs. 

16. The apparatus as set forth in claim 15, wherein the 
means for identifying (384, 500, 504, 508, 514, 530, 534, 538, 
544) vessel center points (402) includes: 

a means for flood-filling (384) imaged vascular structures 
in the slice (376) to form filled regions defined by pixels 
having a first value; and 

a means for iteratively eroding (500, 504, 508, 514, 530, 
534, 538, 544) the edges of the filled regions to identify 
vessel center points (402) . 

17 . The apparatus as set forth in claim 16 further 

including : 

a means for enhancing (78) vessel edge contrast by second 
order spatial differentiation of the image representation (76) 
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to define an edge volume (80) with improved vessel edge 
definition for the flood-filling means (384) to operate upon. 

18. The apparatus as set forth in any one of claims 
16-17, wherein the means for iteratively eroding (500, 504, 
508, 514, 530, 534, 538, 544) the edges of the filled regions 
includes : 

a first erosion pass means for eroding in a first 
direction (500, 504, 508, 514) using a first moving window 
(472) ; 

a second erosion pass means for eroding in a second 
direction (530, 534, 538, 544) using a second moving window 
(474) ; and 

a means for successively invoking the first and second 
erosion pass means (500, 504, 508, 514, 530, 534, 538, 544) a 
plurality of times to iteratively erode the filled region until 
the remaining at least one region (452, 454) corresponds to 

vessel center points. 

19. The apparatus as set forth in any one of claims 
15-18, wherein the means for segmenting, tracking, extracting, 
enhancing, or identifying (630) vascular information (96) 
includes : 

a means for selecting (612) a first vessel center point 
(610) ; 

a means for finding (90) a vessel direction (92) 
corresponding to the first vessel center point (610) based on 
analysis of the image representation (76) in a neighborhood of 
the first vessel center point (610); 

a means for defining a plane (94) in the image 
representation (76) which is perpendicular to the vessel 
direction (92) and contains the first vessel center point 
(610); 
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a means for estimating (86) vessel boundaries (710) 
corresponding to the first vessel center point (610) in the 
defined plane (94) ; 

a means for repeatedly invoking (618, 634) the selecting 
(612), finding (90), defining a plane (94), and estimating (86) 
means for the plurality of vessel center points (84) ; and 

a means for interpolating (770, 774) the estimated vessel 
boundaries (710) to produce a vascular representation. 

20. The apparatus as set forth in claim 19 wherein the 
means for estimating (86) vessel boundaries (710) includes: 

a means for defining (700) an initial geometric contour 
(702) arranged about the vessel center (610) and lying in the 
defined plane (94) ; and 

a means for iteratively optimizing (704) the geometric 
contour constrained to lie in the defined plane (94) . 

21. The apparatus as set forth in either one of claim 19 
and claim 20, wherein the means for finding (90) a vessel 
direction (92) includes: 

a means for constructing (650) a Weingarten matrix (652) ; 

a means for obtaining (654) a plurality of directions 
(656) by eigenvalue decomposition of the Weingarten matrix 
(652) ; and 

a means for selecting (658) the vessel direction (92) from 
the plurality of directions (656) . 

22. The apparatus as set forth in any one of claims 
15-21, further including: 

an intensity inverting means for inverting (156, 160) the 
intensities of the image elements to generate an intensity- 
inverted image (76), the intensity inverting means (156, 160) 
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being invoked conditional (154) upon the angiographic image 
being a black blood angiographic image. 

23 . The apparatus as set forth in any one of claims claim 
15-22, further including: 

a magnetic resonance contrast agent administering means 
(44) for administering a magnetic contrast agent to the patient 
(42) to improve vascular contrast of the three-dimensional 
image representation (72) . 

24 . The apparatus as set forth in any one of claims 
15-23, wherein the imaging means (10) includes at least one of 
a magnetic resonance imaging scanner and a computed tomography 
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